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ABSTRACT 

We simulate a plausible cosmological model in considerable physical and 
numerical detail through the successive phases of reheating (at 10 ^ 2; ^ 20), 
formation of Pop III stars at 2; ~ 15 (due to H2 cooling), with subsequent 
reionization at 2; ~ 7. We assume an efficiency of high mass star formation 
appropriate to leave the universe, after it becomes transparent, with an ionizing 
background J21 ~ 0.4 (at z = 4), near (and perhaps slightly below) the 
observed value. Since the same stars produce the ionizing radiation and the 
first generation of heavy elements, a mean metallicity of (Z/Zq) ~ 1/200 is 
produced in this early phase, but there is a large variation about this mean, 
with the high density regions having Z/Zq ^1/30 and low density regions (or 
the Lyman-alpha forest with A^hi ^ 10^'^'^cm^) having essentially no metals. 

Reionization, when it occurs, is very rapid (phase change-like), which will 
leave a signature which may be detectable by very large area meter-wavelength 
radio instruments. Also, the background UV radiation field will show a sharp 
drop of ~ 10^^ from 1 Ryd to 4 Ryd due to absorption edges. 

The simulated volume is too small to form galaxies, but the smaller 
objects which are found in the simulation obey the Faber- Jackson relation. 

In order to explore theoretically this domain of "the end of the dark ages" 
quantitatively, numerical simulations must have a mass resolution of the order of 
10^-^ Mq in baryons, high spatial resolution (^ 1 kpc to resolve strong clumping, 
and allow for detailed and accurate treatment of both the radiation field and 
atomic/molecular physics. 
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1. Introduction 

In a previous paper (Ostriker & Gnedin 1996; Paper I) we showed, for a specific 
cosmological model, how an early generation of stars would be expected to form from 
the first nonlinear self-gravitating clumps, which developed while the universe was still 
relatively cold, and was primarily filled with neutral gas. This first generation, "Population 
III" , which would reside in widely and fairly smoothly distributed clumps of stars somewhat 
more massive than globular clusters, will reheat and reionize the universe, leading to an 
end of what Martin Rees has called the "Dark Ages" , during the interval 15 > z > 7. Since 
the same stars produce both ionizing photons during their main sequence phase and metals 
when they end their lives as type II supernovae, this early phase of the universe, "the 
clearing of the fog" , also and inevitably leaves an irregular contamination with metals at 
the low (average) level of {Z)/Zq ~ 10"^-^. 

In the current paper we provide the physical basis for the conclusions reached in Paper 
I, and also address a number of related questions, among them being the following: 

1. When, in terms of the Gunn-Peterson effect (Gunn & Peterson 1965), should the 
universe have become transparent? And how closely are the processes labeled 
"reheating" and "reionization" related? 

2. How would one characterize the spatial variations in the metal abundance produced 
by the early generation of stars? 

3. What was the degree of dumpiness in the gas during these early phases? (This is a 
question, typically neglected in most of previous treatments.) 

4. Will specially designed large area telescopes looking at the (rest frame) 21 cm radio 
sky be able to detect this early development of structure? 

5. How important were the various cooling processes (atomic, molecular, and dust), 
the various heating processes (gravitational collapse, UV photons, supernovae), and 
radiative shielding of high density regions (by gas or dust)? 

6. What were the properties of the stellar groups formed at high redshifts? Are they 
observable? 

These and other questions are addressed in the context of a specific cosmological 
model. In a final section we speculate on how robust our conclusions are to variations of 
the assumed model. 
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In discussing various physical effects that play role at the end of the "Dark Ages" , we 
must note the many authors who have studied these processes before, from the pioneering 
work of Couchman & Rees (1986) and Shapiro & Giroux (1987) to the recent careful 
investigations by Shapiro, Giroux, & Babul (1994), Tegmark and Silk (1994, 1995), 
Kawasaki & Fukugita (1994), Shapiro (1995), Giroux & Shapiro (1996), Tegmark et al. 
(1996) to the one-dimensional simulations by Haiman, Thoul, & Loeb (1995) and Haiman 
& Loeb (1996). While varying amounts of physically detailed modelling were included in 
all of those papers, there was a critical factor that could not be followed in any of the 
semi-analytic treatments: the clumping of the gas. However, all the relevant processes are 
dependent on the clumping factor Cbb = Included among these processes are 

recombination, cooling, gravitational collapse, molecular hydrogen formation and numerous 
others. Whereas one-dimensional simulations are able to account for the clumping around 
a single spherically symmetric object, they must neglect the clumping within the flow and 
therefore they cannot predict the average clumping of the gas in the universe, and arc not 
able to follow the evolution of the radiation flcld in any reasonable detail. Thus, even the 
qualitative conclusions of the previous papers may be in doubt until they are confirmed 
by a fully nonlinear three dimensional threatment. A legitimate question to ask at this 
point is the following: Does the present simulation, reported on in this paper have sufficient 
spatial resolution to correctly compute the clumping? The answer is either "yes" or "barely 
enough" , depending on ones degree of optimism. On scales below the Jeans mass clumping 
should be relatively unimportant as it cannot (by dcflnition of the Jeans mass) be assisted 
by gravity. At a characteristic temperature of 10, 000 K, the Jeans mass is around 10^° 
ai z = 4: for the average density of the universe, and is around 10® Mq for the overdensity 
of 10^ (which is reliably resolved by our simulations), whereas our nominal total mass 
resolution is 10^-^ Mq (cf Table 1). Numerical experiments that we have performed indicate 
that the resolution needed to permit the existence of a first (Pop III) generation of stars 
is of the order of 10^'^ Mq, which we just achieve in the computations reported on here. 
Work currently underway should allow us to push the limit to a comoving spatial scale of 
^ 0.5h^^ kpc and a mass scale of 1.0 x 10^ Mq, while simultaneously allowing better for the 
missing large-scale power. 

2. Method 

2.1. Physical Ingredients of a Numerical Simulation 

In order to adequately simulate to at least moderate accuracy the evolution of the 
universe at high redshift, we must allow for all of the most important physical processes. 
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Fortunately, by restricting ourselves to early stages of formation of structure, we can 
incorporate essentially all physics which affects formation of structure at moderate densities 
(n ^ lcm~^), contrary to simulations of galaxy and large-scale structure formation at 
low redshifts, which are generically missing important physical ingredients. The main 
reason behind this difference between simulating high and low rcdshift universe is that 
most of poorly understood complicated physical processes, i.e. effects of heavy elements, 
dust, magnetic fields, are thought to be less important at high redshift, where most of the 
intergalactic gas is of essentially primeval composition and magnetic fields have not yet had 
time to build up to a significant strength (cf Kulsrud et al. 1996). 

We have therefore undertaken to simulate numerically the evolution of the intergalactic 
medium at high redshift, including essentially all physical processes that we identified as 
important for the high redshift evolution of the universe. We have used the SLH-P^M 
cosmological hydrodynamic code as described by Gnedin (1995) and Gnedin & Bertschinger 
(1996) to follow the evolution of the dark matter and the cosmological gas with high 
Lagrangian resolution. We have also included detailed atomic and molecular physics of 
a gas of primeval composition, following in non-equilibrium time-dependent fashion the 
ionization and recombination of hydrogen and helium as well as formation and destruction 
of hydrogen molecules in the ambient radiation field. The radiation field, in turn, allows for 
sources of radiation (quasars and massive stars), sinks (due to continuum opacities) and 
cosmological effects. In particular, we have taken into account the fact that dense lumps will 
be shielded from the background radiation field. This reduces the heating rates for dense 
clumps and makes it nearly certain that once they have formed and started to collapse, 
the process will be irreversible. They will continue to collapse and fragment even though 
outer, lower density regions may reheat and be stabilized against further collapse. However, 
since the precise calculation of the shielding effect requires carrying on the whole radiative 
transfer in the highly inhomogcneous medium, which is well beyond the capabilities of 
existing computers, we adopted a simple approximation which we call the Local Optical 
Depth approximation. In its essentials this approximation treats absorption as localized 
but the sources of radiation as smoothly distributed. We derive radiative transfer equations 
in the expanding universe in the Appendix A, and then introduce the Local Optical Depth 
approximation in the Appendix B. In Appendix C we describe in detail our method of 
following formation and destruction of molecular hydrogen. 

The ionizing radiation is emitted by stars and quasars. However, it is impossible to 
simulate star formation within the cosmological framework directly, and we therefore are 
forced to use a phenomenological approach to account for star formation and feedback 
processes. In regions which are cooling and collapsing we have allowed the formation of 
point-like "stellar" subunits according to the Cen & Ostriker (1992) algorithm as adapted 
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to the SLH code by Gnedin (1996). However, we introduce the following change to the Gen 
& Ostriker algorithm. In their original paper Gen & Ostriker introduced as a numerical 
convenience the density cut-off as a value below which no formation of stellar subunits was 
permitted. The value of the cutoff is arbitrary, and final results are somewhat dependent 
on its value, as has been shown by Gnedin (1996) and Katz, Weinberg, & Hernquist (1996). 
In our improved version of the star formation algorithm, we abandon the density cut-off 
value in the belief that as long as the code has not reached its resolution limit, the evolution 
of the gas is followed accurately, and there is no need to introduce this phenomenology 
altogether. Therefore, we allow for star formation only in the regions which have reached 
the resolution limit, since, in those regions, the further collapse of the gas cannot be followed 
accurately and stellar particles must be created in order to allow the gas to sink beyond the 
resolution limit of a simulation. We identify this limit as the limit when the Lagrangian 
behavior of the SLH code switches to the Eulerian one. Numerically, this corresponds to 
places where the maximum eigenvalue of the softening tensor aij reaches 1/2 (see Gnedin & 
Bertschinger [1996] for the definition of the softening tensor a^), and physically corresponds 
to the spatial resolution shown in Table 1 for the various runs we have performed. These 
numbers are essentially our only free parameters. The sum of them is fixed by matching the 
background radiation field at 2 = 4 and the ratio determines the softness of the ionizing 
background spectrum (cf Fig. |Tl]), also an observable. 

As soon as a stellar particle forms, it releases radiation and (in proportion) metal rich 
gas, which we have considered in the treatment of cooling. Since we cannot distinguish 
between stars and quasars in our treatment of star formation, we allow each stellar particle 
to emit a composite spectrum of radiation, consisting of mixture of both the quasar and 
young star spectra. The shape of the spectrum is shown in Fig. 1 of Gnedin (1996). 
The amount of radiation emitted is again a free parameter of the phenomenological 
star formation algorithm, and we choose it so that we produce ionizing radiation at low 
redshift in the amount which is close to the observed value.0 We assume an efficiency 
^uv,s = 6 X 10^^ for the stellar component of the source spectrum and the equal value of 
^uv,Q = 6 X 10^^ for the QSO-like component of the source spectrum. These numbers are 
essentially our only free parameters. The sum of them is fixed by matching the background 
radiation field at z = 4, and the ratio determined the softness of the ionizing bacground (cf 
Fig. H), also an observable. 

We have also allowed for the injection of the thermal energy from supernova explosions 
into the intergalactic gas according to Gen & Ostriker (1992) methodology. However, at the 



•^However, since we cannot determine the final value for the ionizing intensity until after we have run the 
simulation, we can only approximately reproduce the observational value for the ionizing intensity. 
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resolution we achieve in the simulations reported here the gas density is so high that all this 
thermal energy is immediately radiated away producing negligible net effect on the evolution 
of the universe. This results stems purely from the inadequate treatment of the supernova 
explosions since gas inside the supernova bubbles is too hot to cool. Therefore, the 
appropriate treatment of the supernova explosions would require following hydrodynamics 
of the multi-phase medium (Cowie, McKee, & Ostriker 1981) and modelling interactions 
between the phases, which is beyond the frame of this work. We, therefore, emphasize 
that while supernova explosions are formally incorporated in our simulations, they are not 
followed adequately and, therefore, their effects may be significantly underestimated in our 
treatment. 

The only important piece of physics, that we are totally missing in our simulations, is 
the nonuniformity of the sources of radiation (as explained in Appendix B). One, therefore, 
should bear this in mind when interpreting the results of our simulations. However, since we 
compute all average quantities as appropriate, we expect that the evolution of the universe 
as whole, and, in particular, the evolution of average quantities, is reproduced adequately 
by our simulations. 

2.2. Cosmological Models 

We have adopted a CDM+A cosmological model as a framework for our investigations. 
However, we believe that qualitative results of our simulation are applicable to the whole 
family of the CDM-type models, and we emphasize those qualitative conclusions together 
with specific quantitative results for the model we use. We fix the cosmological parameters 
as follows: 

fio = 0.35, 
h = 0.70, 
Qb = 0.03, 

which is close to the "concordance" model of Ostriker & Steinhardt (1995). We normalize it 
to COBE 10° measurement and allow for a small tilt with n = 0.96. This gives the value of 
(jg = 0.67 for the power spectrum normalization at 8h"^ Mpc at 2; = 0. In order to compute 
the initial conditions accurately, we computed the linear transfer functions for the model 
using the linear Boltzmann code similar to (but different from) the COSMICS package 
(Bertschinger 1995). We show in Fig. |I| transfer functions for both the dark matter (solid 
line) and the baryons (the dotted line) in ratio to the BBKS transfer functions (Bardeen et 
al. 1986) aX z = 100, where we start our simulations. 

We have performed three runs with different box sizes and numerical resolution to 
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Fig. 1. — Transfer functions for the dark matter {solid line) and the baryons {dotted line) aX z = 100 
in ratio to the BBKS transfer function for the CDM+A cosmological model with = 0.35, h = 0.7, 
and = 0.03. 

assess the importance of different scales and estimate the uncertainty due to the finite 
resolution of our simulations. The adopted parameters of those runs are compiled in Table 
1. All runs were stopped ai z — A since at that time the rms density fluctuation at 2/i~^ Mpc 
scales is 0.4 and the absence of nonlinear wavelengths larger than the box size renders the 
simulations senseless. The largest of our simulations, run A, is our fiducial run against 
which to compare results of the two smaller runs. However, it happened by accident that 
the particular realization of initial conditions we used for run A had substantially less power 
on large scales (1 — 2h~^ Mpc) than the true power spectrum, and we have to take this fact 
into account when interpreting our results. In particular, this affects the precise value for 
the redshift of reionization and the evolution at late times, but at earlier times, z ^ 10, 
when the density correlation length is much smaller than the box size, the absence of large 
waves has little effect. 
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Table 1: Numerical Parameters 

Run Box size Total mass res. Spatial res. Dyn. range 

~~^ 128^ 2/1-1 Mpc 3.7 X lO^/i^^M© l.O/i^^kpc 2000 

B 64^ 2/i-iMpc 2.9 X lO^/i-i M0 S.Oh-^kpc 640 

C 64^ 1/i-iMpc 3.7 X lO^/i-^M© l.S/i-^kpc 640 

3. Results 

3.1. Reheating and Reionization 

It is often assumed in the current literature that reheating and reionization of the 
universe at high redshift are two sides of the same physical process, and even the terms 
"reheating" and "reionization" are sometimes used interchangeably. In this section we 
demonstrate that while the physical cause of both of those processes is the same: the high 
energy radiation field, they are nevertheless two different physical processes, separated in 
time as well as proceeding at different speeds. 

Since it is the stars that produce the ionizing radiation, we first show in Figure ^ the 
fraction of the total baryonic mass which is locked in stars (the upper panel) and the 
comoving star formation rate dp^/dt in units of Mq/ yr/ Mpc^, as a function of redshift 
(the middle panel) from out run B. Also in the middle panel we show with the dotted line 
the comoving feedback rate, i.e. the rate of production of high energy photons and metals 
(for convenience, we convert it to the same units). The feedback effects are delayed after 
the formation of a stellar particle according to the methodology of Cen & Ostriker (1992; 
see Gnedin 1996 for details). Note the existence of the first peak of Pop III star formation 
at 2; ~ 16 in otherwise smooth star formation rate. The slow down of the star formation at 
2; < 6 is due to the lack of large-scale waves in our simulation box. 

Fig. ^ shows the temporal evolution for the gas temperature (lower panel) and the 
neutral hydrogen fraction (middle panel). We also show the evolution of the ionizing 
radiation intensity as measured by the quantity J21, called the ionizing intensity, defined as 
follows: 

^ J Jua^du/u ^ 1 

/ Oudu/ V lO^^i erg cm^^ s^^ Hz^^ sr^^ ' 
where is the radiation intensity and Oy is the hydrogen photoionization cross-section. 
The observed values for J21 lie between 0.3 and 1.5 (Savaglio et al. 1996; Cooke, Espey, & 
Carswell 1996) at z ~ 3 - 4. 
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Fig. 2. Upper panel: the fraction os the total baryonic mass locked into stars as a function of 
redshift for run B. Middle panel: the comoving star formation rate in units of M©/ yi/(h~^ Mpc)^ 
{solid line) and the comoving stellar feedback rate, normalized to the same units {dotted line), as 
a function of redshift for the same run. Lower panel: the gas clumping Cbb = {p^)/{p)'^ {solid line) 
and the ionized hydrogen clumping Chu (discused in §3.1; dotted line) as a function for redshift 
for the same run. 

We immediately note that there is a relatively long epoch from ^ = 20 to 2; = 10, 
encompassing many Hubble times at that redshift, when the volume and mass weighted 
average temperature of the universe increases from the value T=llKat2; = 20, dictated 
by the pure adiabatic expansion of the universe, to the value T = 14, 000 K, fixed by the 
atomic physics. This increase in temperature, or reheating of the universe, is accompanied 
by a slow increase in the value of ionizing intensity up to J21 ~ 10~^, but the universe stays 
essentially neutral, with the neutral hydrogen fraction decreasing slowly until z — 7, when 
suddenly the neutral fraction drops by several orders of magnitude at essentially constant 
temperature and the ionizing intensity J21 increases to a value of around 0.5 — 1.0. This 
moment, which happens in a small fraction of a unit redshift (i.e. in a small fraction of 
the then Hubble time), is called reionization of the universe and occurs almost as a phase 
change. 
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Fig. 3. — Evolution of the radiation intensity J21 {upper panel), average neutral hydrogen fraction 
/hi {middle panel), and the average temperature {lower panel) as a function of redshift for run B. 
Solid lines show the volume average and dotted lines show the mass average for the temperature 
and neutral hydrogen fraction. 

We immediately conclude that reheating and reionization are two different physical 
processes; reheating is a slow process happening at high redshift, z > 10, followed by 
a sudden reionization at a lower redshift. During reheating the value of the ionizing 
intensity stays below J21 < 10^'^, because the universe remains mainly neutral at that 
time and the absorption is large. In order to understand why the absorption of ionizing 
photons is significant, let us estimate the recombination time, i.e. the time a hydrogen 
ion needs to recombine after being ionized by a high-energy photon. For the temperature 
T — 10^ K the value of the recombination coefficient R is 4.2 x 10~^^cm^s~^ assuming case 
A recombination. The equation for the ionized hydrogen number density can be written as 
follows: 

^Hii = -2>Hnnn + nniT - Rnenmi, (2) 



- 11 - 



where H is the Hubble constant, F = 4.3 x 10^^^J2is^^ is the photoionization rate, and 
rif, is the free electron number density. Equation (Q) holds at every fluid element in the 
universe. In order to compute the average number density of the ionized hydrogen, we have 
to average equation (0) over the volume of the universe, obtaining: 

?^Hii = -SifnHii + (^Hir) - {Rnerinii), (3) 

where the bar symbol means the volume average. Assuming, that R is constant in space 
and that helium is not ionized (i.e. = riHn), we can obtain the following expression for 
the recombination time: 

^Hii 

or for the ratio of the recombination time to the Hubble time, tn = ^/ H: 



- (4) 



where a is the scale factor, x is the volume average ionization fraction, x = n-nn/nR, and 
Chii is the ionized hydrogen clumping factor, Chu = ('^Hn)/('^Hn)^7 shown in the lower 
panel of Fig. ^ with the dotted line, which is somewhat lower than the total gas clumping 
factor Cbb = because higher density regions are typically less ionized. For the 

model under consideration, h^/^lo = 0.7\/0.35 = 0.41, and for the fully ionized gas, x = 1, 
with no clumping, the recombination time is equal the Hubble time at the redshift Zeq = 8, 
which implies that if the universe stayed uniform, once it was ionized after 2; = 8, it will 
not be able to recombine {t^/tH > 1)- This is a well known result. However, since the 
universe is in fact strongly clumped, and the clumping factor Chu ^ 1, the absorption is 
much higher than that for the uniform universe. In particular, the clumping factor Chu is 
equal 12 at 2; = 8 (and it is even larger in the real universe since we underestimate clumping 
due to the flnite resolution of our simulations), and the recombination time is equal to the 
Hubble time ai z = 8 for the ionization fraction as low as x = I/Chu = 0.08. This means 
that it is sufficient to have only 8% of all hydrogen ionized in order to be able to absorb at 
least one ionizing photon per Hubble time per hydrogen atom (the corresponding number 
is 10% for the recombination case B). 

It is also easy to understand why the reheating phase should precede the reionization 
phase. Let us consider the ratio of the photoheating time, 

^ 3/2kBT 

t-heat — E T ^ 

to the photoionization time, 

_ 1 

^ion = 
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where kB is the Boltzmann constant, and 



Ej = 



is the average energy a hydrogen atom receives in act of ionization. For our case, 
Ej = kB9.3 X 10^ K and T = 10^ K at z = 10, and 

^ = -%^ = 0.16 (6) 

which imphes that the radiation background is much more (by a factor of 6) efficient in 
heating the gas at high redshift than in ionizing it. The reason is simply that, unless the 
ionizing radiation field is very soft, excess heat is liberated for each act of photoionization. 

Finally, we can understand why the reionization is so rapid. Since the recombination 
time is short compared to the Hubble time, there is an approximate balance between the 
emission of ionizing photons from the sources and the absorption (by ionizing a hydrogen 
atom which subsequently recombines). If the emission increases, so does the ionization 
fraction x, the recombination time decreases, and the absorption increases; inversely, the 
decrease in the emission leads to the subsequent decrease in the ionization fraction and an 
increase in the recombination time. Since the star formation rate is increasing with time, as 
shown in the middle panel of Fig. |^, driven by the further collapse of hierarchical structure, 
the ionization fraction slowly increasing in response to the increase in the emission of 
ionizing photons, until it reaches values close to unity. At this point a further increase in 
the recombination rate cannot be achieved, emission and absorption get out of balance, 
and the universe ionizes at the time scale of the emission, which is at the point of loss of 
balance is approximately equal to the recombination time. We, therefore, conclude that the 
sharp decrease in the neutral hydrogen fraction from the values about 0.1 at the beginning 
of (complete) reionization to final values of the order of 10~^ is achieved on a time-scale of 

tH /l + ;zrci\~^/^ 
which for our case of z^ei = 7 has a value of 

. tH 

tre^ ~ TH' 

in excellent agreement with the actual reionization time measured from Fig. ^. This effect of 
fast reionization is entirely due to the high value of the clumping factor. Were the clumping 
ignored in our consideration, we would conclude that the reionization time is more than 10 
times larger than the actual value is, and that reionization is a slow process continuing for 
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more than one Hubble time (as long as the source evolution is also slow as in Fig. This 
is precisely the conclusion achieved in some of the previous work on the reionization where 
the clumping was ignored (see, for example, Haiman & Loeb 1996). 

We point out here that the specifics of the cosmological model does not enter this 
conclusion and the prediction that slow reheating happens first, followed later by sudden 
reionization, is a generic prediction for any cosmological model (with UV photons providing 
both heat and ionization), not necessarily a CDM-type one. 

3.2. The Gunn-Peterson Effect in the Intergalactic Medium 

The complete reionization of the universe manifests itself in the absence of the 
Gunn-Peterson absorption trough (Gunn & Peterson 1965) in the quasar spectra.0 We can 
use the observational limits on the Gunn-Peterson optical depth as a test for cosmological 
models. In order to compute the Gunn-Peterson optical depth, we generate random lines 
of sight through our simulation box|^, and compute the Lyman-alpha absorption r(A) along 
each line of sight. We then compute the average opacity of the intergalactic medium, 

1 ^'"^ 1 f^2^ 

exp(-f ) = ^ E 7,0) (-^(^)) d\ (7) 

where A'ray ~ 1000 is the number of lines of sight, and A^"''' and are wavelength limits 
for the j line of sight. This computation can be compared with observations by Jenkins & 
Ostriker (1991), when expressed as 1 — = exp(— f). 

We show in Figure ^ (the upper panel) the data from Jenkins & Ostriker (1991) as filled 
circles together with our computation (the solid line). We note that our simulation predicts 
somewhat too little absorption compared to the Jenkins & Ostriker data, but the general 
behavior is quite similar to the observations. We reiterate here that we do not expect to 
get an exact match with all existing observational data on high redshift gas contents of the 
universe since we simulate only one cosmological model which may not be correct, and our 
simulations are subject to uncertainties due to the phenomenological description of star 
formation. 



"^Residual clumps (filaments) of neutral gas remain and will produce the Lyman-alpha forest. 

^Those lines of sight are taken to be random with respect to the box orientation, since at lower redshifts, 
when the nonlinear scale approaches the box size, there appear purely artificial correlations between the 
density distribution and the orientation of the computational box. 
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Fig. 4. — Upper panel: continuum suppression factors 1 — Dj^ as a function of rcdshift for Run 
B {solid line) together with the data from Jenkins & Ostriker (1991) (filled circles). Lower panel: 
the Gunn-Peterson optical depth as a function of redshift for Run B (solid line). Shown as a sohd 
triangle is the observational upper limit from Jenkins &; Ostriker (1991) and as an open circle is 
the upper limit from Giallongo et al. (1994). 

We can also express the average absorption along the line of sight as the the average 
Gunn-Peterson optical depth. However, we must exclude the Lyman-alpha forest from the 
total absorption since it is not counted for toward the total Gunn-Peterson optical depth 
(Giallongo et al. 1994). We then compute the average Gunn-Peterson absorption tgp as 
the amount of absorption in the parts of the spectrum that lie above the average opacity, 

1 ^-/a'^' exp(-r(A))^(f-r(A))dA 

exp(-TGp) = J^T. — ' (8) 

j^, e{f-r{X))dX 



•ray j=i 
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where 9{x) is a step function.^ 

We show in the lower panel of Figure ^ the Gunn- Peterson optical depth as a function 
of redshift from our Run B. We also show observational upper limits from Jenkins & 
Ostriker (1991) as a solid triangle and from Giallongo et al. (1994) as an open circle. We 
again note that we predict too little absorption compared to the Jenkins & Ostriker data 
but too much absorption compared to the Giallongo et al. data. Since the Giallongo et al. 
(1994) limit is the most stringent one, we adopt it for our analysis. We note that the result 
of our simulations, while broadly consistent with observations, goes above the observational 
upper limit. This implies that the real universe is more ionized than predicted by our 
simulations. The final value for J21 is 0.34 at z = 4, which is slightly below the currently 
favored observational value as we have noted above. The more realistic value of around 
0.7 — 1.0 would decrease the neutral hydrogen fraction and give lower Gunn-Peterson 
optical depth. However, we cannot rescale the ionizing intensity a posteriori preserving the 
accuracy of our simulations, since the radiation field is tied up in the complex nonlinear 
evolution together with the intergalactic medium^, and we would have to rerun the whole 
simulation with new values for the radiative efficiencies for our stellar particles as explained 
in §2.1. 

Given both the observational and theoretical uncertainties, the agreement between the 
predicted and observed opacity of the intergalactic medium seems quite adequate. It will 
be interesting to see if this rough agreement survives when other cosmological models are 
considered. 



3.3. Metal Enrichment of the Intergalactic Medium 

Since in our simulations we follow the star formation with the phenomenological 
approach, we also account for the dynamics of gas that was processed in stars, and, is, 
therefore, enriched with heavy elements {metals). However, since the IMF at high redshift 
is not known and may deviate significantly from the standard IMF, we have no way to 
assign the value for the metallicity in the processed gas reliably, but we still can study the 



^This method of computing the Gunn-Peterson optical depth is similar to the Miralda-Escude et al. 
(1995) method for identifying Lyman-alpha lines. 

^In particular, even if most of the low density intergalactic medium is in the ionization equilibrium with 
the radiation field, the value of the temperature at every fluid element depends on the previous evolution of 
the density and the radiation field in this fluid element (Miralda-Escude & Rees 1994); since we explicitly 
allow for the nonuniform radiation background, this dependence cannot be parametrized in a simple form. 
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Fig. 5. — The mass weighted {solid line) and the volume weighted {dotted line) fraction of the 
intergalactic gas that was processed in stars as a function of redshift. The dashed hne shows p^/^ 
weighting which roughly corresponds to the weighting by the column density of a Lyman-alpha 
absorber. The right hand scale corresponds to the normalization with the adopted value for the 
yield y = 0.020. 

distribution of this metal enriched gas as compared to the distribution of the primeval gas. 

In Figure ^ we plot the fraction fz of the metal enriched gas in the total intergalactic 
gas as a function of redshift for our fiducial run A. Three different values correspond to 
different weights: the mass weighted value is shown with the solid line, the volume weighted 
value is shown with the dotted line, and the dashed line corresponds to the weighting 
with the p^/^ weight, which roughly reflects the weighting by the column density of a 
Lyman-alpha absorber. Here fz is the ratio of the mass of the gas that was processed in 
stars0 to the total mass of the intergalactic gas. If the "yield" for metal production with 
standard definition is y, than the "metallicity" , or a fraction of mass in elements heavier 
than helium is 9yfz, given our definitions. Thus, for y = 0.020, which is characteristic for 
protoellipticals as measured by the metallicity of the intracluster gas (c.f. Silk 1996), the 
final value of {fz) ~ 5 x 10~^ corresponds to Z 9 x 10~^ or Z/Zq ^ 4.3 x 10~^. For a 
more standard value for the yield {y = 0.005) we would obtain the final metallicity of the 



^We assume here that in average stars loose 10% of their total mass into the intergalactic medimn; this 
amount does not include the mass loss that is confined in the interstellar medium. If this amount is different 
from 10%, the plotted fz should be adjusted simply in proportion to the total fraction of stellar mass lost 
to the intergalactic medium. 
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Fig. 6. — The mass- (upper panel) and volume- {lower panel) weighted probabihty distribution to 
find a fluid element with given values of the fraction of metal enriched gas and the total gas density 
at z = 4. The right hand scale corresponds to the normalization with the adopted value for the 
yield y = 0.020. The bold line shows the average probability, solid contours mark the probability 
above average, and the dashed contours show the probability below average; contour spacing is 
logarithmic with the increment of 1/3. 

order of Z/Zq ^ 1.0 x 10^"^. We also note that hj z = 4 the volume average fraction of the 
metal enriched gas may reach a value of about 0.05%. 

While Fig. ^ shows a rough consistency with observations, it does not show us where 
this metal enriched gas is actually located and it understates the variability in heavy 
element abundances at a given epoch. In Figure ^ we show the probability distribution 
to find a fluid element with the particular values of the fraction of the metal enriched gas 
and the total density of the intergalactic gas, as measured by the overdensity 6 = p/ p — 1. 
The upper panel shows the mass weighted, and the lower panel shows the volume weighted 
distributions. We note that most of the volume has fz well below 0.1%, whereas a 
considerable fraction of mass has fz of the order of 0.1 — 1.0%. We also note that at low 
overdensities, which are characteristic of the low column density Lyman-alpha forest (c.f. 



- 18 - 




Fig. 7. — The average fraction of the metal enriched gas as a function of total gas density at four 
different redshifts z = 4:, 5, 6, and 8. The right hand scale corresponds to the normalization with 
the adopted value for the yield y = 0.020. 

Miralda-Escude et al. 1996; Hui, Gnedin, & Zhang 1996), the range of existing metallicities 
is extremely broad, while the average value is quite low. This is in a qualitative agreement 
with the conclusion by Ranch, Haehnelt, & Steinmetz (1996) who found that observed 
scatter in the abundance ratios of the QSO metal absorption systems exceeds that predicted 
by simulations with the uniform metallicity. 

To quantify this effect further, we plot in Figure ^ the average fraction of the metal 
enriched g function of the total gas density at four different redshifts. In other 

words. Fig. |^ is a Fig. ^ collapsed along the y-direction. It is interesting to note that the 
final values of the metallicity in the high density regions from which galaxies will be made 
approach Z/Zq ^ 1/30, which is close to the minimum metallicity found in the oldest and 
most metal poor Population 11 component of our Galaxy. More that that, this value does 
not change significantly with time after z ~ 8. We also point out that there seems to be 
steady diffusion of heavy elements into low density regions. 

It is remarkable that the average metallicity is only a weak function for wide range of 
overdensities above 6 ~ 10^, and it decreases very fast as the overdensities become smaller. 
This effect allows us to make a specific prediction that the metallicity in the Lyman-alpha 
forest should decrease sharply for the column densities below about iQ^^-^-^"^-^ cm~^ as 
measured ai z = 3 (which roughly corresponds to 5 ~ 30, see Hui et al. 1996). Here we 
allow ourselves a broad error bar in the quoted value of the characteristic column density to 
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account for the possible range of cosmological models and uncertainties in the simulations, 
both numerical uncertainties (like the limitations of the phenomenological star formation 
algorithm), and uncertainties in the input physics (like our lack of knowledge of the IMF at 
high redshift or the precise value of the cosmic baryon density). Recent measurements of 
metal abundances in the Lyman-alpha forest (Songalia & Cowie 1996) found heavy elements 
in a large fraction {'^ 75% of all Lyman-alpha systems with column densities in excess of 
lO^^'^cm"^). Our prediction then implies that, when those studies are extended to lower 
column densities, the sharp decline in the fraction of Lyman-alpha systems with measurable 
heavy elements abundances will be found at column densities of about lO^^'^^^^ *^ cm~^. 
The specific value of the column density at which the metal abundance of the Lyman-alpha 
forest declines can serve as a new test for cosmological models. 

The heavy elements are produced in the regions of star formation, where the gas 
density is high and the potential well is deep, and, therefore one can expect that they are 
effectively trapped inside the (proto-)galaxies, contrary to our finding that they are widely 
distributed even in the low density regions. We must therefore address the question of how 
heavy elements eventually get into the low density regions. Since our simulations do not 
include any special diffusion processes that may be responsible for distributing the metal 
enriched gas over a large range of densities, we conclude that there must exist a mechanism 
which works entirely within the framework of hierarchical clustering. Figure ^ demonstrates 
how this may have happened. We show the distribution of the dark matter (the left 
column), the cosmic gas (the middle column), and the stars (the right column) at three 
different redshifts: z = 5.4 (the upper row), z = 4.9 (the middle row), and z = 4.3 (the 
lower row) represented by particles^ in a thin slice with the width of 12h~^ comoving kpc.^ 
All distances in the slice are in comoving kpc. There are two separate objects at z = 5.4 
marked with the square and the circle that are undergoing merger at lower redshift .0 At 
the lower row aX z = 4.3 the dark matter and the star components of the smaller object 
(circle) passed close to the center of the larger object, while the gas components of both 
objects collided in a gigantic shock and merged at the first impact; as the result the stream 



^While dark matter and stars are numerically represented as particles in the SLH approach, the gas 
component is numerically represented by a distorted quasi-Lagrangian mesh; we convert (for display only) 
cells of this mesh into gas particles for the sake of simplicity and uniform presentation. 

lOWhile dark matter particles all have the same mass, the gas particles and the stellar particles have 
different masses and the visual appearance of the particle density does not necessarily correspond to the 
actual gas or stellar density in the slice. 

^^The objects may have different appearance at different redshifts due to the thinness of the slice: particles 
comprising the objects may enter and leave the slice at different redshifts; we show only a subset of all particles 
in the each object for the sake of clearness. 
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Fig. 8. — Distributions of the dark matter [left column), the cosmic gas {middle column), and the 
stars {right column) at three different redshifts: z = 5.4 {upper row), z = 4.9 {m/iddle row), and 
z = 4.3 {lower row) represented by particles in a thin slice with the width of 12h^^ comoving kpc. 
The square and the circle mark two merging objects, and the oval pointed to by an arrow marks 
ejection of the metal enriched gas from the merger product. 
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of high density gas was ejected from the merged object (marked with the oval in the lower 
middle panel, and pointed to by the arrow). The gas ejected came mostly from the larger 
object, and was heavily enriched by heavy elements. In the next Hubble time or so it will 
be dispersed and mixed with the surrounding primeval intergalactic gas, producing high 
metal abundance even in low density regions lying close to the merged object. But the low 
density gas far from any merger would still be uncontaminated. Thus, the metal enrichment 
of the intergalactic gas will be highly inhomogeneous, as shown in Fig. ^. 

We also note that the fact that we demonstrated a working mechanism for the 
heavy element enrichment does not of course exclude the possibility that there exist other 
mechanisms which distribute metals in the low density intergalactic gas. In particular, 
supernova explosions, which are included in the work of Cen & Ostriker (1992) , are effective 
at dispersing metal enriched gas. 

3.4. The 21cm Emission from High Redshift 

The neutral gas at high redshift is surrounded by the thermal bath of CMB photons, 
which has a temperature much higher than the excitation temperature of the 21 cm line. 
The hydrogen atoms, being excited by the CMB photons, spontaneously decay emitting a 
21 cm line radiation, which being constantly redshifted during the evolution of the universe, 
produces continuum signal at the telescope. As the universe evolves and structure develops, 
velocity focusing in the converging flows will produce enhanced emission at lower redshifts 
until reionization, when the neutral hydrogen fraction drops by several orders of magnitude, 
and 21cm emission declines sharply at the frequency corresponding the 21cm redshifted to 
the epoch of reionization. 

We show in Figure ^ the spectrum of the redshifted 21 cm line as would be observed 
on Earth (the solid line), together with the contribution from the Cosmic Microwave 
background reduced by a factor of a hundred to fit on the same scale (the dotted line). 
The very different spectral signature of the 21 cm radiation makes it possible, at least in 
principle, to distinguish it from the CMB. However, the amplitude of the signal is still 
substantially lower than the sensitivity of all existing or being planned radio telescopes 
(including the Square Kilometer Array). Nevertheless, when (and if) this measurement 
is performed, the reward will be quite substantial: since reionization is so sudden, and, 
therefore, the break in the spectrum at the frequency Ub = ^^0/(1 + -^reion) is also very sharp, 
the measurement of the break frequency will give a very accurate value for the redshift of 
reionization, which, in turn, can place a strong constraint on both cosmological models and 
physical parameters at high redshift. 
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Fig. 9. — The spectrum of the 21 cm emission from the intergalactic gas at high redshift (solid 
line). The sharp break in the spectrum corresponds to the redshift of reionization. The two dashed 
lines show the expected rms fluctuations in this radiation on the 1 {upper curve) and 10 {lower 
curve) arcmin angular scale. Shown with the dotted line is the Cosmic Microwave Background 
radiation intensity reduced by a factor of a hundred to fit on the same scale. 

We also show with the dashed line the rms fluctuations in the 21 cm radiation on 1 and 
10 arcmin angular scale. Since the 1 arcmin corresponds to the comoving scale of around 
1.7h~^ Mpc at z = 7 in the adopted cosmological model, the density fluctuations on this 
scale are of the order of 20%; in the sky those fluctuations will result in the fluctuations 
in the redshifted 21cm radiation. Even if the isotropic component of the 21cm radiation 
is only 1% of the CMB fluctuations, because of the incredible smoothness of the CMB sky, 
AT/T ~ 10~^, fluctuations in the 21cm will dominate the CMB fluctuations by a factor of 
10^ X 10~^ X 0.2 = 200! Since at lower frequencies (earlier times) the density fluctuations 
are smaller, there is a peak in the spectrum of rms fluctuations in the 21 cm radiation close 
to the frequency corresponding to the 21 cm frequency redshifted to z,.ei. This would provide 
an alternative (and, likely, much more efficient) way to measure the redshift of reionization 
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Fig. 10. — The fine structure of tlie 21 cm emission at three different redshifts. Each spectrum is 
centered at the reshifted frequency of the 21 cm. The emission terminates at high velocity values 
due to the final size of the simulation box. 

from the 21 cm from high redshift gas. 

The spectrum in Fig. |^ shows the 21cm emission as measured by a broad filter with 
200 km/ s width. However, since after reionization the most of neutral hydrogen is in 
optically thick lumps (protogalaxies) , the final structure of the 21cm emission at lower 
redshift would consist of a combination of several lines rather than of a continuum emission. 
Figure |l^ serves to illustrate this conclusion. We show the 21 cm emission as a function 
of frequency as measured in velocity units at three different redshifts each centered at the 
frequency corresponding to the redshifted 21cm line. The emission terminates at velocities 
larger than 200 km/ s due to the final size of the simulation box. We note that while at 
high redshift emission is almost uniform in frequency, at lower redshift, after reionization, 
it consists of several overlapping lines. 
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Fig. 11. — The spectra of the background radiation at three different redshifts. Ionization edges 
corresponding to the hydrogen and helium ionization are well observed at all times. 

3.5. Evolution of the Background Radiation 

The radiation emitted by stellar- and quasar-like sources participates in a complex 
interaction with the highly inhomogeneous intergalactic medium. The interplay between 
ionization and recombination processes in the intergalactic gas, coupled with star formation 
and development of cosmic structure determines the shape and the amplitude of the 
background radiation which are not guaranteed to be monotonic in frequency or time. As 
an illustration, we show in Figure |lT] the spectra of the background radiation at three 
different redshifts. This is qualitatively very similar to figure 8 of Ostriker & Cen (1996) 
despite the considerably more detailed treatment of radiative processes in this paper. We 
note that in the 13.6 eV < hu ^ IkeV range the background radiation spectrum has a 
complex shape with pronounced ionization edges at 13.6 eV, 24.6 eV, and 54.4 eV. The 
neutral helium ionization edge at 24.6 eV becomes weak after helium becomes at least singly 
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ionized somewhat prior to hydrogen reionization at z = 7, but the hydrogen ionization 
edge at 13.6 eV is still easily observed even ai z = 4, well after reionization of the universe. 
Examination of Fig. ^ shows that the assumption of a power-law (or other simple form) 
for the background radiation field, which is an approximation adopted by several authors, 
would have been extremely inaccurate due to the strong effects of atomic ionization edges. 

3.6. Properties of Proto-galaxies 

We now concentrate on general properties of objects formed in our simulations, which 
we identify with small proto-galaxies. Paper I addressed in some detail the question of 
the various stellar Populations and star formation rates at different epochs. Here we will 
discuss some of the internal properties of those objects. 

Figure |T2| shows the fraction of metal enriched g clS 3iS Sb function of total mass for 
all objects with the total mass Mtot > W^h-^MQ and stellar mass M* > lO^h ^Mq at 
three different redshifts as marked with different symbols: z = 9.3 (solid squares), z = 5.9 
(stars), and z = 4.3 (open circles). We note that there exists a large fraction of objects with 
as much as 10% of metal enriched gas, which corresponds to solar metallicity if the high 
value of the yield y = 0.020 is adopted as discussed above. We also point out existence of 
an object at 2; ~ 6 with at least the solar metallicity (if the standard value for the yield 
y = 0.005 is adopted) or several times the solar metallicity if the higher value of yield is 
accepted. We therefore conclude that hierarchical clustering models have no difficulty in 
explaining high metallicity observed in quasars at high redshifts. We emphasize here that 
due to the small size of our computational box, our simulations are strongly biased against 
finding very bright massive quasars as those are expected to be rare events in the gaussian 
initial conditions. 

In Figure |13|we show the mass in stars as a fraction of the total mass Mtot (o) and 
the baryonic mass (b) of a given object as a function fo the total mass of the object 
at three different redshifts as described above. A considerable fraction of all objects have 
already turned almost a half of their gas into stars by 2 ~ 4, however, even at the centers 
those objects are still dominated by the dark matter; stars constitute on average one third 
of the total density at the center by 2 ~ 4 for objects with the total mass above 10^/i~^ Mq. 
This is mostly due to our finite resolution {Ih^^ comoving kpc), which corresponds to the 
physical scale of around 300 pc at z ^ 4, and partly due to the fact that objects in Fig. |l^ 
are still efficiently forming stars and actively accreting intergalactic gas which is cooling 
and concentrating at the center of objects, thus increasing the role of baryons compared 
to the dark matter with time, and, therefore, if we were able to continue our simulations 
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Fig. 12. — The fraction of metal enriched gas as a function of the total mass for all bound objects 
with Mtot > W^h-'^MQ and > IO^/i'^Mq at z = 9.3 {solid squares), z = 5.9 (stars), and 
z = 4.3 (open circles). The right hand scale corresponds to the normalization with the adopted 
value for the yield y = 0.020. 

further in time with the appropriate allowance for the missing large-scale waves, we would 
observe the increasing fraction of baryons at the centers of the proto-galaxies. 

Let us now discuss the gravitational structure of proto-galaxies. We show in Figure 
14a the total mass for all objects as a function of dark matter velocity dispersion at three 
different redshifts as explained above. The solid line shows the power-law dependence with 
the power-law index of 2. The fact that our objects follow this power-law shape indicates 



that they are close to the virial equilibrium. In Figure |Tjb we plot the total stellar mass as 
the function of the stellar velocity dispersion, and the solid line is the following law: 



tot 



10^°-^/i~^Mr7 



© 



220 km/ s 



(9) 
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Fig. 13. — The fraction of stars in the total mass (a) and in the total baryonic mass (6) as a function 
of the total mass of an object for all massive objects at three different redshifts as explained before. 

and the power- law index is now 4 instead of 2. If we compare this dependence with the 
TuIIy-Fisher or Faber- Jackson relation we will get a good agreement provided that the 
average stellar mass- to-light ratio is T = = 2.2, which is a reasonable number 

consistent with the solar neighborhood. We also point out that the normalization of our 
"Faber- Jackson" relation is roughly independent of redshift. The scatter is apparently 
substantially larger than the observed scatter in large galaxies (about 1.5-1.7 magnitudes), 
but it mainly arises from the finite mass resolution of our simulations: if the most massive 
objects in our simulations contain many thousands of dark matter particle and even larger 
number of gas particles, they contain a much smaller number of stellar particles.Q We also 
emphasize that the active merging is going on around all massive objects in our simulation 
which tends to increase the scatter. 



10' 



^^The largest object in the simulation has 21,000 dark matter particles, 23,000 gas particles and 2,200 
stellar particles, whereas the objects at the bottom of Fig. |l^ have only 30-40 stellar particles. 
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Fig. 14. — The total mass as a function of the dark matter velocity dispersion (a) and the stellar 
mass as a function of stellar velocity dispersion (&) for all massive objects at three different redshifts 
as explained before. The solid lines have a slope of 2 in (a) and 4 in (b). 

4. Conclusions 

Wc demonstrate by means of a high resolution numerical hydrodynamic simulation how 
reheating and reionization of the universe occur in the CDM+A cosmological model. We 
incorporate in our treatment essentially all relevant physical processes which are important 
for assessing the evolution of average properties of the universe, with one free parameter 
- the efficiency of high mass star formation - determined by matching the computed and 
observed background at redshift z — A. 

We conclude that a generic prediction of models where the universe is reheated 
and reionized by high energy radiation from stellar-like and quasar-like sources is that 
slow reheating happens first followed by sudden reionization. We also demonstrate that 
hierarchical clustering models are broadly consistent (with the appropriate choice of 
parameters) with the measured Gunn-Peterson opacity of the intergalactic medium. 

The average heavy element abundance in the intergalactic medium can reach values 
of 1/200 of solar by 2; = 4 with the large dispersion in the low density regions {6 ^ 30). 
The metal abundances in the high density regions are more uniform and can reach 1/30 
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in average, but there exist individual bound object (proto-galaxies) that can have solar 
or higher abundance of heavy elements at z ~ 6. This fact implies that hierarchical 
clustering models like the one studied in this paper can easily explain observed high 
metallicity of quasars at z = 5. We also propose a working mechanism based on mergers for 
transporting the heavy elements, initially locked in the high density regions, into the lower 
density intergalactic medium, where they are observed in (sufficiently high column density) 
Lyman-alpha absorbers. 

Reionization, which suddenly reduces the neutral hydrogen abundance by several 
orders of magnitude, leaves a distinct signature in the rcdshiftcd 21cm emission from the 
intergalactic hydrogen. While this emission cannot be measured by existing or proposed 
meter-wavelength radio telescopes, which lack sensitivity by a factor of several hundred, 
future ultra-high sensitivity experiments may detect this emission, giving us a very strong 
(due to the suddenness of reionization) constraint on the evolution of the universe at high 
redshift. 

The Faber-Jackson relation, the correlation between the luminosity and the velocity 
dispersion of a galaxy is reproduced by our simulations (assuming constant mass-to-light 
ratio) with correct amplitude and slope, indicating that this phenomenon can be explained 
entirely based on hierarchical clustering model and simple cooling arguments. Both, the 
slope and the amplitude of our Faber-Jackson relation, are independent of redshift within 
numerical uncertainties. 
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A. Radiative Transfer Equations in the Expanding Universe 

Let us start with deriving the full radiative transfer equation in comoving coordinates. 
Let F{t,x\p'') be the distribution function for photons in comoving coordinates and 
comoving momentum 

k k 
c 

where a = a{t) is the scale factor, h is the Planck constant, u is the photon frequency, and 
n'^ is the unit vector in the direction of photon propagation. By definition (Peebles 1980), 



N^^ J F{t,x\p'')d^xd^p 
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is the number of photons in the universe; A^^ changes only due to photon absorption or 
emission in physical processes (ionization, recombination, cooling, etc). We can, therefore, 
write down the continuity equation for the function F in the phase-space {x\p^): 

where the term RHS(F) denotes photon absorption and emission in physical processes. 

It is customary, however, to use specific intensity Jy{t^x\n^) of radiation instead of 
F{t,x^,p'^), where, by definition, J^dv dVLdAdt is the energy of photons with frequences 
from V to V + dv passing through the area dA in the time interval dt in the solid angle dVL 
around . It is easy to convert from F to recalling that d^x = cdtdA and d^p = p^dpdfl: 

d^xd^p 3^^i^^ 
du dQ dA dt (? 
Substituting F from (|A2|) into ([Al|), we easily obtain: 



J, = hvF ^"Z./. . =a'-^F (A2) 



2 



- 3- + 77- + — ^-.TT- P'p-TTT-,Ju = RHS(J.), (A3) 



dt a dx^ ^ ' p^ dp \ a^h^u^ 

where we have written the divergence over momentum in spherical coordinates explicitly 
presenting the derivative with respect to the total momentum and noting that the 
derivatives with respect to orientation vanish due to the isotropicity of the Robertson- 
Walker space-time. Since p = ahu/c, we finally obtain the following most general radiative 
transfer equation: 

where H is the Hubble constant, H = d/a, and we specified physical processes including 
absorption with the coefficient ki, and emission (sources) Si/. 

We now consider the mean specific intensity in the universe, averaged over all space 
and all directions, 

lit) = {Mt,x\n^)), (A5) 

where the averaging operator () acting on a function /(x^n'^) of position and direction is 

defined as: ^ 

{f(x\n^)) = a number = lim / d^x / dVLfix\n''). (A6) 

v^oo AttV Jv J 

Applying the averaging operator to the equation ( |A4|) , we obtain the following equation for 
the mean specific intensity: 

^-H (u^ - 3 J.) = -hJ. + S., (A7) 
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where, by definition, 
and 

= = . 

We would like to stress here that we use the bar symbol merely as a part of notation 
and not to denote the space average; in particular, is not a space average of k^ since it is 
averaged, weighted by the local value of specific intensity J^. 



The purpose of this section is, however, not to derive equation (|A7|) , but an equation 
for fluctuations in the specific intensity around the mean. It is convenient therefore to 
introduce a new quantity called relative (to the mean) specific intensity, fy{t,x\n^), as a 
measure of deviation of the specific intensity in a particular position in space x* and in a 
particular direction n'^: 

Ju = fuJy, 

SO that {fy) = 1. It is straightforward to derive the following equation for fy. 



dt (9x* ^ ' du 

where 



+ {i'fu) = - {ky - ky)fy + y , (A9) 



= ^^V^. (MO) 



Equation (|A9|) is tedious and by no means simpler than the original equation ([A4|) . 
However, if we restrict ourselves to scales significantly smaller than the horizon size, and 
matter velocities much smaller than the speed of light (Newtonian limit), then we can 
simplify equation ([A9| ) substantially by noting that the relative specific intensity does not 
change substantially and the universe does not expand substantially over the period of time 
a photon needs to travel over the scale under consideration. Therefore, we can ignore the 
time derivative and the derivative with respect to the frequency. Under these assumptions, 
equation ( |A9| ) reduces to the following simple equation: 

c df 

-^'iri = -{ky~ky)fy + ijy (All) 

a ax* 

Equations ([A7| ) and ( |A11| ) are represent the radiative transfer equations in the expanding 
universe in the Newtonian limit. 
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B. Local Optical Depth Approximation to the Radiative Transfer Equations 

Our goal in this section is to describe numerical techniques for solving equations ( [A7| ) 
and (JAIID within the framework of simulating large scale structure and galaxy formation in 



the universe. The former of those two equations can be reduced to an ordinary differential 
equation and offers no difficulties in solving it (cf Gnedin 1996). The latter one however 
is beyond practical implementation for the current computer capabilities and resolution 
requirements for galaxy formation simulations. 

The main difficulty comes from the fact that fi, is a function of six variables: three 
spatial coordinates and three coordinates in the momentum space. In addition, unless 
Ji, by itself is needed in a simulation, only direction-averaged quantities are required to 
calculate photoionization and photoheating rates which are needed to accurately compute 
gas temperature, ionization fractions, cooling rates etc. But even if we approximate the 
directional dependence of fi,, it still will be of little use since f^^ will remain a function of 
four arguments; storing such a function with sufficient number of resolution elements in 
each of four direction is a challenging problem for modern galaxy formation simulations. 

The second complication with practical implementation of the radiative transfer 
equation ( |A11 ) is that the ionization state of cosmic gas (and, therefore, the absorption 



coefficient k^) can change on a very short time-scale, substantially shorter than the 
hydrodynamic Courant time-scale. Equation (|A11|) must therefore be solved millions of 



times for a typical large cosmological simulation, which is obviously well beyond existing 
computer power. 

Fortunately, there exist two physical reasons why a reasonable approximation to 
equation ( |A11D can be found, that will employ only functions of three arguments. 

There are two main physical effects which produce spatially variable specific intensity: 
the optical depth of the gas, that can shield high density regions from external radiations 
{ku is a function of position), and nonuniformly distributed sources of radiation (■j/'j, is 
a function of position). The optical depth is indeed a function of seven arguments: the 
frequency z/ and two spatial positions (optical depth from one point to another); however, its 
dependence on frequency is uniquely determined since for known column densities of neutral 
hydrogen A^hi, neutral helium N^ei, and singly ionized helium Nn^u (we consider only 
hydrogen and helium contributions to the optical depth since heavy element abundances 
are small) the optical depth is 

Tu = Nu,a^' + iVHci^?^^ + Nucu(T^''\ (Bl) 
where a* is the photoionization cross-section for a species i = H I, He I, He ll and is a known 
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function of frequency. In addition, there is no need to know with high accuracy, since if 
the optical depth is small, it does not matter if it is 10"'^ or 10^^, and if it is large, it again 
does not matter if it is 10^ or 10^. It is only r,^ ~ 1 that should be computed accurately. 
However, since ~ 1 primarily happens at the boundary of an optically thick region, the 
total volume occupied by regions with ~ 1 is not large and little error in physical state 
of gas is introduced from a large error in the optical depth calculation |^. 

The second important physical reason why equation (|A11| ) can be reasonably 
approximated with functions of positions only, is that source function ipi, dependence on 
frequency, if known, does not change with time since photons does not have enough time to 
redshift significantly while traveling a distance short compared to the horizon scale. Since 
ipi, includes only fluctuations in the source function, it is mainly contributed to by a few 
nearby sources, while the redshift of radiation emitted by distant (and, therefore, uniform) 
sources is already taken into account by a derivative with respect to v in equation (|^). 
Thus, the frequency dependence of relative specific intensity at every particular point 
and at every moment in time can be easily calculated from scratch and needs not to be 
saved as an additional dimension for fy (in particular, it is the same for all positions). 

Let us now proceed further to deriving the simplest possible approximation that has 
any sense, which we call the Local Optical Depth Approximation. First, let us consider the 
homogeneous equation (|A11|) when no spatial dependence of sources is taken into account. 



Let f^^ be its solution. It is convenient and instructive to introduce optical depth t^, as: 

/f(t,x\n'=) = exp (- [Ty{t,x\n^) - f,(t)]) (B2) 
where fy{t) is defined as: 



because {fy) = 1. Substituting ( p2D into (|A11|) with ijjy = Q V4e obtain the following 
equation: 

n'^^^=-%-K) (B3) 
with the formal solution (omitting the time dependency for clearness): 

Ty{x\n^) = r^°)(x' - nXxV),n^) + - T \K{x' + nH) - K] dl, (B4) 

where Tj^"^ is an arbitrary function of two two-dimensional vectors which is specified from 
boundary conditions. One can see that is indeed an optical depth in a sense that it is an 
integral of an absorption coefficient over a distance. 



^•^Those handwaving arguments can not serve as a proof to the approximation described below but merely 
as justification for searching for reasonable approximations. 
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However, as we discussed before, there is no practical sense to implement this equation 
in a real code, or even in a direction- averaged version of the exact solution (in fact, it is 
not even clear how to calculate the direction-averaged /jy in a closed form from ||B4|]). More 



than that, what we would eventually like to get as an approximation to fi^ is the completely 
local dependence of on properties of the gas; this would allow us to consider cooling 
evolution within a hydrodynamic time-step of each resolution element independently and 
use an ordinary differential equation solver to resolve smallest time-scales in the problem 
(see Gnedin 1996 for details). We therefore introduce an approximation to fi, by defining 
Tj, in (|B2|) as a function of position and frequency only by the following ansatz: 

T,{x') = [nu,a^' + nneia?^' + nHen^?^"] L, (B5) 

where nm, nnei, and nneii are local number densities of neutral hydrogen, neutral helium, 
and singly ionized helium and L is a characteristic distance (which is a function of position, 
but not direction; the resulting fjf^ is therefore a function of position and frequency, 
but not direction). The justification to this ansatz is that it is likely that ionization 
fractions are similar in the neighboring points and if we choose the characteristic distance L 
appropriately, we will obtain a successful approximation to f^. Since L is obviously related 
to the characteristic scale over which the density changes, we choose the following form for 
L: 

P (B6) 



a|Vp|V/3p|Ap| 



where parameters a = 2.11 and (3 = 0.27 are fitted to produce the correct value for the 
column density for a spherical density distribution with inverse square law profile in two 
limits r — > and r ^ oo (two constraints produce two parameter values); these values for a 
and P also give a very close match to the gaussian density profile and for r""^ (photoionized 
inverse square law) profiles. 

The Local Optical Depth approximation is uncontrolled approximation in a sense that 
it seems impossible to define its range of applicability; we can guarantee that it is not zero 
since it gives correct results for a spherical density distribution with the inverse square law 
profile and quite accurate results for other density profiles. In addition, we have tested the 
approximation by solving equation ( |A11| ) numerically for density distributions extracted at 



a few fixed epochs from simulations discussed in this paper; the correspondence between 
the exact result and the approximate /j, is remarkably good. The epochs, however, have 
been chosen after the universe reionizes, so, fu deviates from unity appreciably only within 
isolated dense optically thick clouds and the Local Optical Depth approximation seems to 
give reasonable results in this limit; when the most of the universe is "optically thick" (i.e. 
when > 1) the approximation is less accurate, but since r,^ > 1 in most of the volume 
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of the universe in that case, big errors in t^, do not lead to substantial errors in ionization 
abundances as have been argued above. 

Now, equipped with the approximation to the solution /^^^ for the homogeneous 
radiative transfer equation, we can obtain solution for the whole equation ( |A11|) . Since the 
direction enters equation ( [All|) only as a parameter, we can introduce a distance I along 
the direction and equation ( |A11| ) becomes an ordinary differential equation, which can 
be easily solved (we assume further that ip^ does not depend on direction, i.e. that radiation 
sources are isotropic; in the opposite case dependence on direction can not be eliminated 
and the described approximation becomes invalid): 



c Jo 

where we specify t„{x\x\) below. Since no directional dependence required in order to 
calculate effects of radiation on the gas, we average equation (p7| ) over all directions; in 
order to keep notation compact, we retain the same symbol fu{x^) for the direction averaged 
f^{x\n^), i.e. 

f^(x')^j-IUx\n')dn. 

Then equation ( P7| ) becomes: 



If we introduce a new vector x\ = + n^l and notice that / = |a;* — x^l and d^xi = f'dl dQ, 
we obtain the following equation for fy{x^): 



c J (x* — x\)^ 
where Tu{x\x\) is the true optical depth from point x* to point x\. 

Integral ( p9|) , except for the exponential factor, is a potential for the force with 
force law, which can be solved for by known numerical techniques, for example, P^M. The 
exponential factor, however, offers a considerably difficulty in solving integral (|B9D faster 
than in 0{N'^) operations. We, therefore, leave this form of the solution to the radiative 
transfer equation ( |A11| ) until further work. 



C. Molecular Hydrogen 



In this section we describe our method for solving for nonequilibrium time-dependent 
molecular hydrogen abundance. The method we adopt is exact but assumes that molecular 
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hydrogen abundance is small compared to the neutral hydrogen abundance and, thus, 
we can follow the ionization balance without taking molecular hydrogen formation into 
account. 

We use Shapiro & Kang (1987) paper as our source of reaction rates. The Table | 
contains all relevant reactions together with their number from Shapiro & Kang (1987) 
paper; we use this number as an identification for a reaction. 

Assuming that we know time-dependence of hydrogen and helium ion abundances as a 
function of time, we can write down the following three equations for three chemical species 
under consideration: H2, Hg", and H~: 



1 dXn- 
Uh dt 



-Xh- [kgXm + heXe + hrXm + hsX^n + hgXnn + k2iX^+ + P27) + 
(A;7Xh:X, + A;i3Xh,X,), (Cla) 



1 dX^+ 

' - {kwXm + ^20^6 + hiXn- + P28 + Pso) + 



and 



rib dt 

(/cqXhiXhii + /C12XH2XH11 + A;i9XhiiXh- + P29XH2) , (Clb) 



1 (iXn 

^ = -Xh2 (A;iiXhi + A;i2XHn + hzXe + h^Xe + fclS^Ha + P29 + P3l) + 



Tib dt 

(/csXhiXh- + A;ioXhiXh+ + A;2iXh+Xh-) , (Clc) 
where the fractional abundance Xj for a species j is defined as: 

=X. — , (C2) 

where itih is the mass of a hydrogen atom and p is the gas mass density; the total baryon 
number density is 

= (C3) 
and temperature and Xj for j = H I, H ll. He I, He ll. He ill, e are known functions of time. 



Equations ^I] have disadvantage that they are nonlinear, and, therefore, do not offer 
a simple solution. The nonlinearity comes in with two terms: ki^X^^ and A;2iXjj+Xh-. 
However, there is a lucky coincidence that each of those terms is superseded by another 
term in each of equations (|C1|) : 



1. fc2iXjj+XH- ^ /cisXhiiXh- in equation for H since ^21 ~ kis; 
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2. A;2iXjj+Xh- <^ k2oX^+Xe in equation for H2 since ^21 ~ ^20! 

3. k2iX^+Xii- ^ ksXuiXii- and ki^X^^ <^ /chXhj^hi in equation for H2 since 
^21 ^ 10'^ /es and /cis ~ /cn, provided that X^+ -C 10~^Xhi (this condition is always 
satisfied in simulations described in this paper). 



Neglecting two nonlinear terms, we get a system of three linear ordinary differential 
equations which offers an analytical solution. 

Let us introduce a three-dimensional vector Y as follows: 



Y 



(C4) 



V Xn, J 

Then equations ( pi| ) can be written in a matrix form as follows: 

Y = -AY + S, 

where A is a matrix, 



(C5) 



A = nb 



f Di + k^Xui + A^igXHii 

— kigX-^u 

—ksXni 



-A;i3Xe \ 

D2 + kioX-^i —kuX-iiu — P29 

-kioXui + ki2Xuu + ki^Xe + P29 / 



(C6) 



D is a vector, 



D 



+ P2I \ 

k2oXe + P28 + P30 

^ kiiXm + kuXe + P31 J 



(C7) 



and 5 is a vector, 



S = Tib 



( kjXn.X, \ 
kgXu iXh 11 




(a 



Equation ( |U5D can be solved analytically in a time interval At from to to ti = to + At 
in which both A and S do not change significantly with the O(At^) to obtain: 



where 



(C9) 



(CIO) 
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and 



At 



to 



S{t)dt. 



(Cll) 



Equation (CDf) then gives the molecular hydrogen and and H fractions at the time 
moment ti. 
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Table 2: Molecular Hydrogen Chemistry Rates 

No.*^ Reaction Rate notation 

7 H I + e ^ H- k'j 

8 H I + H- H2 + e ks 

9 H I + H II ^ kg 

10 Hl + H+^Ha + Hll k^o 

11 Hl + H2^3Hl ku 

12 H2 + H II ^ H I + ki2 

13 H2 + e ^ H I + H- A;i3 

14 H2 + e ^ 2H I + e ku 

15 H2 + H2^H2 + 2Hl fci5 

16 H- + e ^ H I + 2e fcig 

17 H- + Hi^2Hl + e A;i7 

18 H- + Hll^2Hl kis 

19 H- + Hii^H^ + e A;i9 

20 + e ^ 2H I k2o 

21 H^ + H-^H2 + Hl k2i 

27 H- + 7 ^ H I + e A;27 

28 + 7 ^ H I + H II k28 

29 H2 + 7 H+ + e fc29 

30 + 7 ^ 2H II + e k^o 

31 H2 + 7^2Hl fcsi 



"as in Sliapiro & Kang 1987 



